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ABSTRACT 

A  new  algorithm  for  identifying  and  characterizing  vortices  in  complex  flows  is  presented. 
The  scheme  uses  both  the  vorticity  and  pressure  fields.  A  skeleton  line  along  the  center  of  a 
vortex  is  produced  by  a  two-step  predictor-corrector  scheme.  The  technique  uses  the  vector 
field  to  move  in  the  direction  of  the  skeleton  line  and  the  scalar  field  to  correct  the  location  in 
the  plane  perpendicular  to  the  skeleton  line.  A  general  vortex  cross  section  can  be  concisely 
defined  with  five  parameters  at  each  point  along  the  skeleton  line.  The  details  of  the  method 
and  examples  of  its  use  are  discussed. 
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1  Introduction 


Vortices  are  considered  the  most  important  structures  that  control  tlie  dynamics  of  flow 
flelds.  Large-scale  vortices  observed  in  atmospheric  and  oceanographic  flows  are  responsible 
for  ozone  holes,  hurricanes,  tornadoes,  and  maelstroms.  Vortices  that  are  shed  from  aircraft 
canards,  wings,  and  control  surfaces  affect  the  handling  characteristics  of  the  airplane  and 
determine  the  closest  safe  following  distance  of  other  aircraft.  On  a  smaller  scale,  vortices  are 
the  fundamental  building  blocks  of  turbulent  flow  structures  and  are  intimately  connected 
with  the  mixing  of  chemical  species,  heat  transfer,  and  drag  forces. 

Although  the  term  “vortex”  connotes  a  similar  concept  in  the  minds  of  most  fluid  dynam- 
icists,  a  precise  definition  is  difficult  to  obtain.  Robinson  [1]  suggests  the  following  working 
definition:  “A  vortex  exists  when  instantaneous  streamlines  mapped  onto  a  plane  normal  to 
the  vortex  core  exhibit  a  roughly  circular  or  spiral  patlerti,  when  viewed  from  a  reference 
frame  moving  with  the  center  of  the  vortex  core.”  Unfortunately,  this  definition  requires  a 
knowledge  of  the  vortex  core  before  one  can  determine  whether  something  is  a  vortex.  In 
practice,  Robinson  (2]  and  Robinson,  Kline,  and  Spalart  [3]  use  the  above  rigorous  definition 
to  confirm  that  a  particular  structure  is,  in  fact,  a  vortex.  Regions  of  low  pressure  are  used 
to  identify  candidate  vortices.  Their  scheme  exploits  the  fact  that  the  pressure  in  the  core  of 
a  vortex  is  lower  than  that  of  the  surrounding  fluid.  A  radial  force  is  needed  to  provide  the 
centripetal  acceleration  that  keeps  a  particle  rotating  about  an  axis.  In  the  case  of  vortices, 
this  force  is  in  the  form  of  the  pressure;  the  pressure  inside  a  vortex  is  lower  than  the  pressure 
outside  the  vortex.  Robinson  and  his  colleagues  find  that  elongated  low-pressure  regions  in 
incompressible  turbulent  flows  almost  always  indicate  vortex  cores,  but  these  surfaces  can 
be  difficult  to  characterize  and  provide  no  information  in  regard  to  the  sense  of  rotation  or 
the  connectivity  between  the  structures. 

Moin  and  Kim  [4]  and  Kim  and  Moin  [5]  use  vorticity  lines  (sometimes  called  vorticity 
field  lines  or  vortex  lines)  to  visualize  vortical  structures  in  turbulent  channel  flow.  A  vorticity 
line  is  everywhere  parallel  to  the  local  vorticity  vector  and  is  defined  by 
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where  tD  =  V  x  u  is  the  vorticity  vector,  u  is  the  velocity  vector,  x  is  the  location  in  space,  and 
s  is  the  distance  along  the  vorticity  line.  The  definition  suggests  that  vorticity  lines  should 
go  through  the  cores  of  vortices,  which  is  usually  the  Ccise.  However,  the  resulting  vortex 
lines  are  extremely  sensitive  to  the  choice  of  initial  location  for  the  integration.  As  Moin  and 
Kim  [4]  point  out,  “If  we  choose  Xq  arbitrarily,  the  resulting  vortex  line  is  likely  to  wander 
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over  the  whole  flow  field  like  a  badly  tangled  fishing  line,  and  it  would  be  very  difficult  to  ^ 
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identify  the  organized  structures  (if  any)  through  which  the  line  may  have  passed.”  Kim  and 
Moin  [5]  illustrate  the  potential  tangle  in  their  Fig.  2.  To  avoid  such  a  confusing  jumble,  they 
carefully  select  the  initial  points  (inside  suspected  vortices)  and  focus  on  specific  parts  of 
the  resultant  lines  to  illustrate  some  important  features  of  turbulent  channel  flow.  However, 
Robinson  [1]  shows  that  even  experienced  researchers  can  be  surprisingly  misled  by  ordinary 
vorticity  lines. 

Villasenor  and  Vincent  [6]  present  an  algorithm  for  the  recognition  and  visualization  of 
elongated  vortices  in  three-dimensional  time-dependent  flow  fields.  They  start  from  a  seed 
point  and  compute  the  average  length  of  all  vorticity  vectors  contained  in  a  small-radius 
cylinder.  The  axis  of  the  cylinder  extends  from  the  seed  point  to  the  surface  of  a  sphere 
centered  about  the  seed  point.  They  repeat  this  step  for  a  large  number  of  cylinders  that 
emanate  from  the  seed  point.  The  cylinder  with  the  maximum  average  vorticity  vector  length 
becomes  a  segment  of  the  vortex  tube.  This  process  is  repeated  until  the  vorticity  decrecises 
enough  to  end  the  vortex  tube.  Because  the  scheme  depends  only  on  the  magnitude  of 
the  vorticity  vectors  rather  than  their  direction,  the  algorithm  must  limit  the  angle  between 
adjacent  line  segments  so  that  the  new  line  segment  is  different  from  the  previous  line  segment 
traced  in  the  opposite  direction.  All  vorticity  vectors  inside  the  vortex  tube  are  drawn.  The 
process  can  be  repeated  for  multiple  time  steps  so  that  the  evolution  of  the  vortex  tubes  is 
visualized  as  a  sequence  of  all  stored  vorticity  vectors,  Villasenor  and  Vincent  [6]  claim  that 
the  use  of  only  the  magnitudes  and  not  the  directions  of  the  vorticity  vectors  is  advantageous 
because  it  allows  the  algorithm  to  find  structures  that  do  not  have  axes  aligned  with  the 
principal  axis  of  the  structure.  As  a  consequence,  the  algorithm  can  inadvertently  capture 
structures  that  are  not  vortices.  In  addition,  by  visualizing  all  vectors  that  are  inside  the 
cylinder,  structures  are  included  that  are  not  part  of  the  vortex  originally  intended  to  be 
viewed.  An  example  can  be  seen  in  Fig.  4(c)  of  Villasenor  and  Vincent  [6],  in  which  another 
structure  is  responsible  for  some  of  the  vorticity  vectors  that  are  almost  orthogonal  to  the 
vortex  axis. 

Zabusky,  Boratav,  Pelz,  Gao,  Silver,  and  Cooper  [7]  fit  ellipsoids  to  surfaces  of  constant 
vorticity  |u;|  and  to  constant  vortex  stretching  |u;  •  Vu|/(u;|  in  an  effort  to  understand  the 
dynamics  of  a  vortex  reconnection  process.  A  bundle  of  vector  field  lines  of  the  vorticity 
u5  and  the  vortex  stretching  cD  •  Vu  emanate  from  the  major  and  minor  axes  and  from  the 
center  of  the  minor  ellipse  of  the  respective  fitted  ellipsoids.  Both  the  ellipsoids  and  the 
vector  field  lines  provide  useful  information  about  the  vortices  for  their  flow  field,  which 
contains  neither  solid  boundaries  nor  a  mean  straining  field.  In  flows  with  these  additional 
complications,  the  regions  with  large  vorticity  magnitudes  do  not  necessarily  correspond  to 
vortices;  hence,  ellipsoids  of  constant  vorticity  coupled  with  vorticity  lines  are  unlikely  to 
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provide  useful  characterizations. 

Chong,  Perry,  and  Cantwell  [8]  address  the  question  of  when  a  region  of  vorticity  is 
a  vortex.  They  suggest  that  a  vortex  core  is  a  region  where  the  velocity-gradient  tensor 
has  complex  eigenvalues.  In  such  a  region,  the  rotation  tensor  dominates  over  the  rate-of- 
strain  tensor.  Soria  and  Cantwell  [9]  use  this  approach  to  study  vortical  structures  in  free- 
shear  flows.  At  those  points  at  which  the  vorticity  magnitude  exceeds  a  specified  fraction 
of  the  maximum  vorticity,  the  eigenvalues  of  the  velocity-gradient  matrix  are  determined. 
Those  regions  that  have  complex  eigenvalues  and  satisfy  the  vorticity-threshold  criterion  arc 
rendered  as  solid  surfaces  that  represent  the  vortices  in  the  flow.  This  method  correctly 
identifies  the  large  vortical  structures  in  the  flow.  However,  the  method  clutters  the  picture 
with  many  smaller  structures  and  does  not  provide  a  way  to  link  the  smaller  vortical  volumes 
with  the  larger  coherent  vortices  of  which  they  might  be  a  part. 

Yates  and  Chapman  [10]  carefully  explore  two  definitions  of  vortex  cores.  One  definition 
associates  a  vortex  core  with  a  streamline  of  minimum  curvature  within  a  region  of  spiraling 
streamlines.  The  other  definition  regards  the  vortex  core  as  the  line  defined  by  the  local 
maxima  of  normalized  helicity.  Under  certain  circumstances,  both  definitions  produce  the 
same  vortex  core.  Unfortunately,  all  of  the  analyses  and  conclusions  are  appropriate  only 
for  steady  flows. 

Bernard,  Thomas,  and  Handler  [11]  use  a  semiautomated  procedure  to  identify  quasi- 
streamwise  vortices  in  the  near-wall  region  of  turbulent-channel-flow  direct-numerical-simulation 
data.  Their  method  finds  local  centers  of  rotation  in  user-specified  regions  in  planes  perpen¬ 
dicular  to  the  streamwise  direction.  The  local  centers  of  rotation  are  linked  to  corresponding 
centers  in  adjacent  streamwise-normal  planes  to  form  a  vortex-core  line.  The  results  of  their 
work  indicate  that  experienced  users  can  correctly  find  essentially  all  critical  vortices  re¬ 
sponsible  for  the  maintenance  of  the  Reynolds  stress.  However,  their  method  fails  to  capture 
vortices  with  axes  that  are  not  nearly  aligned  with  the  streamwise  direction.  This  shortcom¬ 
ing  can  be  a  major  difficulty  for  flows  such  as  free-shear  layers  and  transitional  boundary 
layers,  in  which  both  streamwise-  and  spanwise-oriented  vortices  are  important.  In  addition, 
the  procedure  depends  heavily  on  user  intervention  to  indicate  regions  where  the  program 
should  search  for  centers  of  rotation.  This  process  is  tedious,  and  the  detailed  results  depend 
upon  nonquantifiable  whims  of  the  user. 

Here,  we  present  a  new  approach  for  identifying  and  characterizing  vortices  in  complex 
flow  fields.  Rather  than  a  dependence  upon  a  single  quantity  to  define  the  vortex,  the 
new  approach  is  a  predictor-corrector  method  that  converges  to  the  vortex  skeletons  (the 
lines  that  pass  through  the  vortex  cores).  In  section  2,  we  describe  how  the  new  method 
determines  the  vortex  skeletons.  In  section  3,  we  discuss  some  novel  vortex  cross-section 
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schemes  that  put  meat  on  the  skeleton  lines,  and  we  summarize  our  findings  in  section  4. 

2  The  vorticity-predictor  pressure-corrector  method 

Our  new  predictor-corrector  method  produces  an  ordered  set  of  points  that  approximates 
the  skeleton  of  a  vortex.  Associated  with  each  point  are  quantities  that  describe*  the  local 
characteristics  of  the  vortex.  These  quantities  may  include  the  vorticity,  the  pressure,  cross- 
sectional  information,  the  circulation,  and  other  quantities  of  interest. 

The  new  method  produces  lines  that  are  similar  to  vorticity  lines.  However,  because  a 
distinction  exists  between  the  vorticity,  which  is  a  mathematical  function  of  the  instantaneous 
velocity  field,  and  a  vortex,  which  is  a  physical  structure  with  coherence  over  a  region  of 
space,  the  straightforward  integration  of  the  vorticity  lines  must  be  modified  to  locate  the 
skeleton  lines  associated  with  individual  vortices.  At  the  very  least,  because  vorticity  lines 
can  begin  and  end  only  at  domain  boundaries  and  vortices  have  no  such  restriction,  we  must 
develop  criteria  for  initiating  and  terminating  the  integration.  In  addition,  the  method  must 
be  self-correcting  (i.e.,  line  trajectories  that  diverge  from  the  vortex  core  should  be  subject 
to  a  restoring  (correcting)  force). 

In  this  section,  we  will  first  discuss  the  procedure  used  to  find  an  initial  point  on  the 
vortex-core  skeleton.  Then  we  will  reveal  the  details  involved  in  growing  the  vortex-core 
skeleton  from  the  seed  point.  Finally,  the  termination  of  the  vortex  skeleton  will  be  addressed. 

2.1  Finding  a  seed  point 

Any  flow-field  information  can  be  used  to  find  starting  locations  for  the  predictor-corrector 
process.  The  pressure  and  the  vorticity  magnitude  are  convenient  scalars  for  choosing  seed 
points.  In  our  work,  the  flow  field  is  scanned  in  planes  perpendicular  to  the  streamwise 
direction.  The  scanning  direction  affects  the  order  in  which  vortices  are  located,  but  not  the 
overall  features  of  the  vortices.  In  each  plane,  the  values  of  the  pressure  and  the  vorticity 
magnitude  are  checked.  Threshold  values  can  be  chosen  a  priori,  or  they  can  be  a  predeter¬ 
mined  fraction  of  the  extrema.  If  the  pressure  is  less  than  the  pressure  threshold  and  the 
vorticity  magnitude  is  greater  than  its  threshold,  then  the  point  location  is  further  refined. 
In  our  implementation,  this  positional  refinement  shifts  the  seed  point  to  the  location  of  the 
local  pressure  minimum  in  the  plane  perpendicular  to  the  vorticity  vector.  (The  rationale 
for  this  refinement  will  be  discussed  later.)  The  sampling  can  be  done  at  every  grid  point. 
We  experimented  with  sampling  at  every  second,  third,  and  fourth  points.  Although  a  test 
of  each  grid  point  to  determine  if  the  specified  pressure  and  vorticity  criteria  are  met  is  not 
expensive,  we  had  hoped  that  sampling  for  seed  points  on  a  coarser  mesh  would  eliminate 
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the  multitude  of  redundant  vortices  that  we  found.  Global  coarsening  of  the  grid  does  not 
uniformly  eliminate  the  problem;  some  of  the  redundancies  are  avoided,  but  others  remain. 

A  better  method  for  eliminating  the  redundant  vortices  involves  the  selective  elimination 
of  grid  points  from  the  supply  of  possible  seeds.  Elements  of  a  three-dimensional  integer 
array  correspond  to  volumes  about  each  grid  point.  All  of  the  elements  are  originally  set 
to  zero.  For  the  first  vortex,  any  point  in  the  flow  field  may  be  evaluated  as  a  potential 
seed  point.  As  points  are  added  to  the  vortex  skeleton,  the  elements  of  the  integer  array 
associated  with  the  points  are  (lagged.  The  flagged  points  will  not  be  available  eis  seeds 
for  subsequent  vortices.  Although  many  coincident  vortices  are  eliminated  in  this  way,  a 
further  reduction  can  be  achieved  by  also  flagging  elements  of  the  array  that  correspond  to 
points  that  are  within  the  vortex  cross  section.  In  Fig.  1,  multiple  skeleton  lines  started  from 
different  seed  points  all  pass  through  the  same  vortex  tube.  The  redundancies  are  eliminated 
when  points  inside  a  tube  are  excluded  from  the  pool  of  future  seed  points. 

2.2  Growing  the  skeleton 

Once  a  seed  point  has  been  selected,  the  skeleton  of  the  vortex  core  can  be  grown  from 
this  point.  The  scheme  that  we  have  developed  is  a  two-stage  predictor-corrector  method. 
With  this  technique,  the  next  position  of  the  vortex  skeleton  is  predicted  by  integrating 
along  the  vorticity  vector.  This  candidate  location  is  corrected  by  adjusting  the  position 
to  the  pressure  minimum  in  the  plane  that  is  perpendicular  to  the  vorticity  vector.  To 
ensure  that  the  minimum  is  actually  part  of  the  vortex  under  consideration,  only  a  limited 
angle  between  the  vorticity  vectors  at  the  predicted  and  the  corrected  point  is  allowed.  ’ 
The  continuous  modification  of  the  vorticity-line  location  lessens  the  sensitivity  to  both  the 
initial  conditions  and  the  integration  details.  These  sensitivities  are  common  problems  in 
vorticity-line  calculations.  Just  as  importantly,  the  corrector  stage  decreases  the  influence 
of  background  mean  shear  and  small-scale  vorticity  fluctuations  on  the  identification  of  the 
vortex  cores.  The  identified  vortex  skeleton  corresponds  more  closely  to  what  one  might 
intuitively  choose  as  the  vortex  skeleton  than  to  the  standard  vorticity-line  result. 

The  modified  vortex-line  scheme  is  illustrated  in  the  schematic  diagrams  of  Fig.  2  and  is 
also  summarized  in  pseudocode  below. 


'This  differs  from  the  angle  limitation  imposed  by  Villasenor  and  Vincent  (6).  They  restrict  the  change 
in  direction  of  the  skeleton  line  along  the  center  of  the  vortex  tube  from  one  point  on  the  skeleton  line  to 
the  next.  We  restrict  the  change  in  the  direction  of  the  vorticity  vector  from  the  predicted  location  to  the 
corrected  location. 


For  each  remaining  seed  point 


If  the  seed  point  po  is  not  in  any  previous  vortex 

While  the  vortex  skeleton  continues 

1.  Determine  vorticity  cji  at  p,- 

2.  Integrate  vorticity  to  find  p^+i  (predictor) 

3.  Determine  vorticity  aJj+i  at  pj+i  (corrector) 

4.  Find  location  p,>i  of  minimum  pressure  in 
plane  perpendicular  to  at  pi+i 

5.  Calculate  quantities  of  interest  at  p{+i 

6 .  t  i —  t  "f"  1 

The  calculation  of  the  vortex  skeleton  proceeds  both  forward  and  backward  from  the  seed 
point.  The  details  for  continuing  the  calculation  from  one  point  to  the  next  are  indicated 
by  the  numbered  items  in  the  pseudocode.  Processes  1  and  2,  which  are  shown  in  Fig.  2(a), 
represent  the  predictor  stage  of  the  algorithm.  The  corrector  stage  is  summarized  by  pro¬ 
cesses  3  and  4,  which  are  illustrated  in  Figs.  2(b)  and  2(c).  We  iterate  the  corrector  stage  so 
that  the  pressure  at  position  pi+i  is  the  local  minimum  in  the  plane  perpendicular  to  Wi+i. 
We  restrict  the  movement  in  the  corrector  steps  by  limiting  the  angle  between  and  a5,+i. 
The  resultant  state,  which  is  illustrated  in  Fig.  2(d),  is  equivalent  to  that  in  Fig.  2(a),  except 
that  the  index  *  -I- 1  replaces  the  index  i. 

The  effectiveness  of  the  predictor-corrector  scheme  is  illustrated  in  Fig.  3,  in  which  data 
from  the  direct  numerical  simulations  of  Singer  and  Joslin  [12]  are  analyzed.  The  transparent 
vortex  tube  is  constructed  with  data  from  the  full  predictor-corrector  method.  Its  core  is 
indicated  by  the  darker  skeleton.  The  lighter  skeleton  follows  the  uncorrected  integral  curve 
of  the  vorticity.  It  is  obtained  by  disabling  the  corrector  phase  of  the  scheme.  The  vorticity 
line  deviates  from  the  core,  exits  the  vortex  tube  entirely,  and  wanders  in  the  flow  field. 

Although  the  general  behavior  of  the  predictor-corrector  algorithm  is  reliable  and  robust, 
optimal  performance  of  the  technique  requires  careful  attention  to  implementation  details. 
The  remainder  of  this  subsection  addresses  issues  that  are  important  to  the  successful  use 
of  this  method. 

2.2.1  Eliminating  feeder  vortices 

Because  the  vorticity  near  the  edge  of  a  vortex  may  be  skewed  with  respect  to  the  vorticity  at 
the  vortex  center,  the  location  of  the  pressure  minimum  in  the  plane  perpendicular  to  the  edge 
vorticity  might  not  be  in  the  vortex  center.  This  potential  mismatch  of  the  pressure-minimum 


and  the  vortex-center  locations  is  rarely  a  problem  in  the  vortex  skeleton-line  calculations 
because  the  predictor-corrector  method  inhibits  large  excursions  from  the  vortex  center. 
However,  a  seed  point  might  be  selected  near  the  vortex  edge.  In  this  case,  adjustment  of 
the  point  location  to  the  pressure  minimum  might  still  result  in  a  point  that  is  far  from  the 
vortex  center.  These  situations  result  in  small  “feeder  vortices”  that  spiral  toward  the  vortex 
center.  Examples  of  feeder  vortices  are  illustrated  in  Fig.  4.  We  found  that  we  can  eliminate 
most  of  these  feeder  vortices  by  taking  advantage  of  the  asymmetry  of  the  predictor-corrector 
method.  In  either  forward  or  backward  mode,  the  predictor-corrector  scheme  will  converge 
to  a  vortex  center.  We  can  eliminate  the  feeder  vortices  by  redefining  the  initial  point  (i.e., 
the  point  at  which  we  start  recording  the  elements  on  the  skeleton  line)  to  be  the  location 
after  a  fixed  number  (usually  5  to  10)  of  predictor-corrector  steps.  If  the  original  seed  point 
is  located  on  the  vortex  core,  the  method  will  return  to  the  seed  on  the  backward  integration 
path.  However,  if  the  original  seed  point  is  located  on  a  feeder  vortex,  then  the  backward 
integration  will  miss  the  original  seed  point  and  the  skeleton  of  the  feeder  vortex.  Instead, 
the  retraced  skeleton  line  will  closely  follow  the  vortex  center.  Although  the  use  of  many 
integration  points  before  the  initial  point  will  eliminate  feeder  vortices,  this  practice  will  also 
limit  the  minimum  extent  of  any  detected  vortices. 

2.2.2  Numerical  considerations 

Neither  the  pressure  minimum  nor  the  result  of  the  vorticity-line  integration  is  likely  to 
be  on  a  grid  point;  hence,  we  must  choose  an  interpolation  scheme  to  find  the  pressure 
and  vorticity  at  arbitrary  locations  in  the  flow  field.  Our  first  experiments  with  trilinear 
interpolation  were  a  disappointment.  No  local  extrema  could  be  found  between  two  points 
connected  by  a  straight  line.  To  allow  the  local  pressure  minima  to  exist  in  the  interior  of  a 
grid  cell,  a  higher  order  interpolation  method  is  necessary.  Second-order  interpolation  would 
use  three  grid  points;  two  from  one  side  of  the  desired  location  and  one  from  the  other  side. 
To  reduce  any  biais  from  the  interpolation,  a  four-point  Lagrange  interpolation  is  used  in  each 
of  the  three  coordinate  directions.  When  the  outer  points  of  the  four-point  method  are  not 
in  the  domain,  two-point  Lagrange  (linear)  interpolation  is  used.  The  complete  interpolation 
scheme  works  quite  well,  although  it  requires  more  computer  time  than  any  other  subroutine 
in  the  computer  program. 

The  vorticity  integration  can  be  performed  with  a  variety  of  methods.  First,  we  used  a 
fourth-order  Runge-Kutta  approach.  This  produced  satisfactory  results;  however,  step-size 
optimization  was  difficult  to  automate.  Instead,  we  developed  a  new  technique  whereby  the 
point-to-point  distance  in  the  vorticity  integration  is  always  equal  to  the  smallest  dimension 
of  the  local  grid  cell.  The  new  point  location  is  found  by  advancing  this  distance  in  the 
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direction  of  the  local  vorticity  vector.  This  procedure  ensures  that  successive  points  will  not 
be  more  than  one  grid  cell  apart,  so  that  if  the  original  calculation  is  well  resolved,  then  the 
vorticity- line  calculation  will  also  be  sufficiently  resolved.  The  procedure  also  reduces  the 
chance  of  wasting  many  calculations  inside  a  single  grid  cell. 

2.2.3  Pressure-minimum  corrector  step 

The  pressure-minimum  correction  scheme  uses  the  method  of  steepest  descent  to  find  the 
local  pressure  minimum  in  the  plane  perpendicular  to  the  vorticity  vector.  First,  the  vorticity 
vector  at  the  candidate  point  is  determined  by  interpolation  from  the  surrounding  grid  points. 
Two  perpendicular  unit  vectors  /  and  g  (Fig.  5)  are  determined  in  the  plane  normal  to  the 
vorticity  vector.  The  smallest  grid-cell  dimension  d  is  used  as  a  local  length  scale  to  find 

A 

two  points  in  the  plane  in  the  /  and  g  directions.  The  pressures  at  these  two  points  and 
the  current  point  are  used  to  form  the  in-plane  pressure-gradient  vector.  A  sample  point 
is  chosen  at  a  distance  d  from  the  current  point  in  the  direction  opposite  to  the  pressure- 
gradient  vector.  If  the  pressure  at  the  new  sample  point  is  less  than  the  pressure  at  the 
current  point  and  the  dot  product  of  the  vorticities  at  the  sample  and  current  points  is 
greater  than  a  specified  value  (0.90-0.95  works  well),  then  the  new  sample  point  becomes 
the  current  location  and  the  process  is  repeated.  If  the  sample  point  heis  a  greater  pressure 
than  the  current  point,  then  the  distance  d  is  halved  and  a  new  sample-point  comparison 
is  made.  If  the  orginal  point  is  still  the  point  of  minimum  pressure  after  d  has  been  halved 
twice,  then  one  last  check  for  a  possible  minimum  is  made  in  the  direction  of  the  positive 
gradient.  (Rare  circumstances  with  a  great  deal  of  symmetry  can  produce  a  pressure-gradient 
vector  that  is  in  the  wrong  direction.)  Failure  to  find  any  points  with  lower  pressure  leads 
to  acceptance  of  the  current  point  as  the  next  point  on  the  vortex-core  skeleton. 

2.3  The  end  of  the  line 

Vorticity  lines  extend  until  they  intersect  a  domain  boundary,  but  vortices  typically  begin 
and  end  inside  the  domain.  Here,  we  discuss  some  techniques  for  stopping  the  integra¬ 
tion.  One  particularly  clean  termination  occurs  when  the  vortex  cross  section,  which  will 
be  discussed  in  the  next  section,  has  zero  area.  This  approach  provides  smooth  surfaces  for 
visualization  with  no  abrupt  vortex-tube  cutoITs  in  the  direction  that  is  tangential  to  the 
vorticity.  Although  we  have  used  this  line-termination  technique  successfully,  the  method 
can  fail  to  show  connections  between  parts  of  the  same  vortex.  For  instance,  if  a  low-intensity 
region  exists  between  high-intensity  regions  of  the  same  vortex,  then  the  low-intensity  region 
might  not  satisfy  the  criteria  for  a  finite  cross  section.  If  both  high-intensity  regions  have 


8 


finite  cross  sections,  then  the  single  vortex  can  be  educed  as  two  small  disconnected  vortices. 
Although  the  criteria  that  define  the  lateral  edge  of  the  vortex  can  be  weakened  to  ensure 
that  the  low-intensity  region  has  a  finite  cross  section,  the  problem  can  potentially  reemerge 
with  a  new  vortex  that  is  educed  with  the  weaker  cross-section  criteria.  Our  resolution  of 
this  problem  exploits  the  asymmetric  nature  of  the  predictor-correction  method. 

Because  the  predictor-corrector  method  will  follow  the  core  of  a  vortex  regardless  of  the 
criteria  used  to  define  the  vortex  cross  section,  the  vortex  skeleton  line  is  continued  even 
when  the  cross-sectional  area  of  the  vortex  is  equal  to  zero.  The  vortex  of  interest  can 
either  reintensify  or  dissipate.  If  the  vortex  reintensifies,  then  the  continuation  of  the  vortex 
skeleton  line  will  provide  a  link  between  the  two  more  intense  regions  of  the  vortex.  This 
link  can  be  visualized  as  a  thread  that  connects  the  two  disjoint  regions,  or  the  two  regions 
can  simply  be  rendered  with  the  same  color  or  texture.  On  the  other  hand,  if  the  vortex 
dissipates,  then  the  vortex  skeleton  line  continuation  will  wander  through  the  flow  field  and 
eventually  either  intercept  a  domain  boundary  or  enter  a  new  vortex.  If  a  domain  boundary 
is  reached,  then  the  elements  of  the  vortex  skeleton  line  that  were  computed  after  the  cross- 
sectional  area  became  zero  are  discarded.  These  same  points  are  also  discarded  if  the  vortex 
skeleton  line  enters  a  new  vortex  with  a  nonzero  cross-sectional  area.  To  determine  whether 
the  new  region  of  finite  cross  section  is  a  continuation  of  the  original  vortex  or  an  entirely  new 
vortex,  we  march  the  predictor-corrector  scheme  backwards  for  the  same  number  of  steps 
taken  since  the  previous  region  of  nonzero  cross  section  was  exited.  Some  possible  scenarios 
are  illustrated  in  Fig.  6.  In  Fig.  6(a),  the  skeleton  line  leaves  the  first  vortex  tube  at  point 
Pi  and  continues  for  N  steps  until  it  encounters  the  second  vortex  tube  at  point  ^2-  The 
predictor-corrector  scheme  is  then  marched  backwards  N  steps  from  P2  to  P3.  The  distance 
between  points  Pi  and  P3  is  small  relative  to  the  distance  between  Pi  and  P2  (a  10-percent 
criterion  is  used);  hence,  the  link  between  Pi  and  P2  is  most  probably  a  low-intensity  vortex, 
and  we  keep  the  thread  between  these  vortex  tubes.  However,  in  Fig.  6(b)  the  vortex  tube 
that  ends  at  P4  continues  to  dissipate,  and  the  continuation  of  its  skeleton  line  lacks  clear 
direction  and  wanders  through  the  flow  field.  The  line  intercepts  another  vortex  tube  at  P5 
after  M  steps.  The  predictor-corrector  method  is  marched  backwards  M  steps  from  P5  to 
Pq.  Initially,  the  reverse  integration  retraces  the  forward  integration,  but  halfway  between 
Ps  and  Pe  the  two  lines  diverge  rapidly  and  become  uncorrelated.  The  distance  from  P4 
to  Pe  is  a  large  fraction  of  the  distance  from  P4  to  P5,  so  the  algorithm  concludes  that  the 
vortex  tube  intersected  at  P5  is  different  from  the  vortex  tube  that  ends  at  P4.  The  points 
on  the  vortex  skeleton  line  that  connect  the  two  tubes  are  discarded,  and  the  vortex  skeleton 
is  terminated.  Finally,  in  Fig.  6(c),  the  continuation  of  the  skeleton  line  of  the  vortex  tube 
that  ends  at  point  P?  intersects  the  side  of  another  vortex  tube  (shown  as  a  wireframe)  and 
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is  immediately  taken  to  the  pressure  minimum  at  fg-  The  reverse  integration  for  this  case 
is  along  the  axis  of  the  new  vortex  tube  away  from  the  original  vortex.  The  point  Pg  is  far 
from  Pr;  hence,  the  two  vortex  tubes  are  distinct  from  each  other  and  the  line  connecting 
them  is  discarded. 

3  Putting  meat  on  the  bones 

The  determination  of  an  appropriate  vortex  cross  section  hcis  been  one  of  the  more  difficult 
practical  tispects  of  this  work.  For  isolated  vortices,  a  simple  pressure  criterion  to  define  the 
edge  of  a  vortex  works  quite  well,  although  the  information  content  of  the  visualization  is 
little  more  than  that  which  is  available  from  pressure  isosurfaces.  In  regions  where  vortices 
interact,  the  pressure  alone  is  inadequate.  The  low-pressure  regions  from  two  or  more  vortices 
can  merge  and  distort  the  radius  estimate  of  any  single  vortex.  A  similar  problem  arises  if  the 
vortex  edge  is  defined  in  terms  of  vorticity  magnitude.  This  particular  difficulty  is  resolved 
if  the  angle  between  the  vorticity  vector  on  the  skeleton  line  and  the  vorticity  vector  at  any 
radial  position  is  restricted.  Any  angle  greater  than  90®  indicates  that  the  fluid  at  the  radial 
position  is  rotating  in  the  direction  opposite  to  that  in  the  core.  The  vortex  circulation  F  is 
defined  as 

r  =  .  kdA  (2) 

where  n  is  a  unit  vector  normal  to  the  surface  S]  dA  is  an  area  element  on  the  surface 
S\  and  5  is  a  cross-sectional  surface  of  the  vortex.  The  vortex  circulation  increases  as  the 
cross-sectional  area  increases,  provided  that  the  angle  between  the  vorticity  vector  at  a  given 
radial  position  and  the  vorticity  vector  on  the  skeleton  line  is  less  than  90®.  This  observation 
suggests  the  use  of  the  90®  variation  in  vorticity  vectors  as  a  single  vortex-edge  criterion. 
Unfortunately,  in  a  uniform  shear  flow,  this  criterion  results  in  an  infinite  cross  section 
because  the  vorticity  at  all  radial  locations  is  the  same.  In  practice,  we  have  found  that  the 
90®  restriction  works  well  in  combination  with  a  low-pressure  criterion  for  the  vortex  edge. 
For  the  actual  computation  of  the  radial  distance,  the  pressure  and  the  vorticity  must  be 
sampled  along  a  number  of  radial  lines  perpendicular  to  the  vortex  skeleton  line.  We  sample 
at  an  increment  that  is  equal  to  the  minimum  grid-cell  dimension  at  the  skeleton  line.  The 
dot  product  of  the  vorticity  at  the  radial  point  with  the  vorticity  on  the  skeleton  line  must  be 
greater  than  zero  for  the  90®-variation  criterion  to  hold.  At  each  position,  this  dot  product 
and  the  pressure  are  checked  to  ensure  that  they  both  satisfy  the  respective  requirements.  If 
either  fails,  the  radial  position  along  that  line  is  the  point  of  failure  as  determined  by  linear 
interpolation  between  the  two  most  recent  sample  points.  By  sampling  along  radial  lines 
that  emanate  from  a  point  on  the  skeleton  line,  we  implicitly  require  that  the  cross  section 
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be  starlike  (i.e.,  each  radial  line  will  intersect  the  cross-section  boundary  only  once.)  Cross 
sections  with  more  complicated  geometries  are  truncated  to  a  starlike  form. 

Relative  measures  for  the  vortex  edge  are  also  possible.  A  relative  measure  defines  the 
vortex  edge  as  the  position  at  which  some  quantity  decreases  below  a  specified  fraction  of  its 
value  on  the  skeleton  line.  Clearly,  this  measure  will  always  give  a  finite  cross  section,  and 
the  cross  section  can  become  infinite  in  a  uniform  shear  flow.  We  have  found  that  relative 
measures  are  not  as  informative  as  the  absolute  measures. 

I'he  characterization  of  the  cross  section  can  take  many  forms.  The  detailed  aspects  of 
the  data  are  most  nearly  recovered  if  the  radial  location  of  the  vortex  edge  is  retained  at 
many  closely  spaced  azimuthal  angles.  If  the  radius  of  the  cross  section  were  sampled  at 
10°  increments  in  the  azimuthal  direction,  then  36  radial  distances  and  a  reference  vector 
at  every  position  along  the  vortex  skeleton  line  would  need  to  be  stored.  For  the  minimal- 
storage  approach  that  provides  characteristic  cross-sectional  information,  all  of  the  radial 
distances  are  averaged  to  obtain  a  single  scalar  radius.  No  reference  vector  is  needed  if  only 
the  radius  is  used,  so  this  simple  method  uses  the  least  amount  of  storage.  We  have  found 
that  the  radius  is  a  very  good  cross-section  descriptor  for  isolated  vortices.  When  vortices 
begin  to  interact,  the  radius  does  not  provide  a  good  description  of  the  asymmetric  distortion 
experienced  by  the  vortices.  However,  the  first  few  coefficients  of  the  Fourier  series  of  the 
radial  locations  provides  a  convenient  compromise  between  the  simple  radius  approach  and 
a  full  description  of  all  radial  locations.  The  series  are  easy  to  compute,  easy  to  interpret, 
and  allow  a  large  range  of  cross-sectional  shapes.  In  our  work,  we  keep  the  constant  term 
(which  is  the  radius),  the  first  and  second  sine  and  cosine  coefficients,  and  a  reference  vector. 
Most  of  the  cases  that  we  have  checked  have  a  factor-of-10  drop  in  the  magnitude  of  the 
first  and  second  coefficients,  a  fact  that  suggests  that  the  series  is  convergent.  Figure  7 
illustrates  a  single  cross  section  of  a  vortex  educed  from  direct  numerical  simulation  data. 
The  shaded  region  represents  the  interior  of  the  vortex  tube,  the  boundary  of  which  was 
determined  by  sampling  along  lines  that  radiate  from  the  cross  (which  marks  the  skeleton¬ 
line  location)  at  1°  intervals.  The.  dashed  line  is  a  circle  centered  about  the  skeleton-line 
location  with  the  averaged  radius  of  the  vortex  tube.  The  solid  line  is  the  two-term  Fourier 
series  representation  of  the  vortex  cross  section.  Note  how  the  vortex  boundary  represented 
by  the  Fourier  series  more  nearly  shows  both  the  eccentricity  and  the  flattening  of  the  vortex 
tube  compared  to  the  circular  cross  section. 
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4  Conclusions 


The  innovative  use  of  a  two-step  predictor-corrector  algorithm  heis  been  introduced  to  educe 
vortices  from  flow-field  data.  Rather  than  relying  on  any  single  quantity  to  determine  the 
vortex  skeleton  line,  the  new  method  uses  vorticity  to  predict  the  new  location  of  the  next 
element  along  the  line  and  then  corrects  this  location  by  using  the  pressure  field  in  the 
plane  perpendicular  to  the  vorticity  vector.  The  prediction  stage  resembles  a  portion  of 
a  vorticity-line  calculation;  the  correction  stage  maintains  the  line  near  the  vortex  center. 
Unlike  other  approaches,  our  method  is  able  to  treat  the  skeleton  line  through  the  vortex 
core  as  an  attractor  in  the  flow  field. 

To  make  the  method  work  in  all  but  the  simplest  flow  fifelds,  numerous  side  issues, 
which  are  common  to  all  techniques  that  grow  a  skeleton  line  from  seed  points,  must  be 
addressed.  This  paper  discusses  a  number  of  novel  approaches  that  we  have  developed  to 
deal  with  matters  such  as  seed-point  selection,  feeder-vortex  elimination,  vortex  skeleton¬ 
line  termination,  and  vortex  cross-section  description.  Sample  extractions  of  vortices  from 
various  flow  fields  illustrate  the  different  aspects  of  the  technique. 

The  predictor-corrector  technique  presented  here  can  identify  vortices  in  flow  fields  that 
are  far  more  complex  than  those  that  we  have  used  here  for  illustration  purposes.  An 
extension  of  the  direct  numerical  simulation  of  Singer  and  Joslin  [12]  provides  one  of  the  more 
interesting  flow  configurations  that  we  have  considered.  The  predictor-corrector  algorithm 
educes  a  complex  tangle  of  many  interconnected  vortices  that  are  visualized  in  Fig.  8.  All  of 
the  implementation  features  discussed  above  are  used  to  identify  the  vortices  in  this  messy 
flow  field. 

Although  many  modifications  can  be  made  to  the  basic  scheme,  the  principle  of  using  a 
vector  field  to  predict  the  location  of  the  next  point  and  a  scalar  field  to  correct  this  position 
distinguishes  this  method  from  others. 
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Figure  1:  Multiple  realizations  of  same  vortex  tube  from  different  seed  points. 


Figure  2:  Schematic  of  predictor-corrector  algorithm,  (a)  Vorticity  tDj  is  used  to  predict 
candidate  point  p,+i.  (b)  Pressure  is  determined  in  plane  perpendicular  to  vorticity  uJi+i. 
(c)  New  point  on  skeleton  line  p,+i  is  corrected  to  pressure  minimum  in  plane  perpendicular 
to  vorticity  Wi+i.  (d)  Repeat  at  point  pi+i. 
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Figure  3:  Standard  vorticity  lines  are  compared  with  results  from  predictor-corrector  scheme. 
Both  results  use  same  flow  data  and  same  computer  code.  Standard  vorticity  lines  are 
obtained  by  disabling  corrector  portion  of  program.  Transparent  tube  is  constructed  from 
data  from  predictor-corrector  method. 


Figure  4:  Feeder  vortices  merge  with  a  large-scale  hairpin  vortex.  Points  that  satisfy  seed- 
point  criteria  exist  on  edge  of  vortex  tube.  They  curve  inward  toward  centerline  and  then 
follow  main  skeleton  line. 
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Figure  5:  Position  correction  using  local  pressure  field.  Orthogonal  unit  vectors  /  and  g  are 
in  plane  perpendicular  to  vorticity  vector.  Method  of  steepest  descent  is  used  to  locate  local 
pressure  minimum. 
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Figure  6(a). 


Figure  6(c). 


Figure  6:  Schematic  of  line-termination  procedure,  (a)  Forward  integration  from  P\  to  Pi 
gives  approximately  same  path  as  reverse  integration  from  Pi  to  P3,  so  points  Pi  and  Pi 
are  connected  by  weak  vortex,  (b)  Forward  integration  from  P4  to  Ps  differs  markedly  from 
reverse  integration  from  Ps  to  Ps,  so  vortex  terminates  at  P4.  (c)  Forward  integration  from 
Pj  to  Pg  intersects  the  side  of  new  vortex  tube  (shown  as  wireframe).  Reverse  integration 
from  Pa  to  Pg  goes  along  axis  of  new  vortex  away  from  original  tube.  Distance  from  P7  to 
Pg  is  large,  so  two  vortices  are  not  linked. 
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Figure  7:  Comparison  of  different  ways  to  represent  vortex  cross  section.  Central  cross  in¬ 
dicates  point  on  skeleton  line.  Finely  sampled  vortex  cross  section  is  represented  by  shaded 
region.  Dashed  line  is  circle  with  averaged  radius.  Solid  line  is  two-term  Fourier  representa¬ 
tion  of  vortex  cross  section. 
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predictor-corrector  algorithm. 
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